insurance <- read.csv("insurance_cost.csv")
summary(insurance)
## age sex bmi children
## Min. :18.00 Length:1338 Min. :15.96 Min. :0.000
## 1st Qu.:27.00 Class :character 1st Qu.:26.30 1st Qu.:0.000
## Median :39.00 Mode :character Median :30.40 Median :1.000
## Mean :39.21 Mean :30.66 Mean :1.095
## 3rd Qu.:51.00 3rd Qu.:34.69 3rd Qu.:2.000
## Max. :64.00 Max. :53.13 Max. :5.000
## smoker region charges
## Length:1338 Length:1338 Min. : 1122
## Class :character Class :character 1st Qu.: 4740
## Mode :character Mode :character Median : 9382
## Mean :13270
## 3rd Qu.:16640
## Max. :63770
library(dplyr)
library(ggplot2)
library(plotly)
library(ggbiplot)
library(ggpubr)
library(corrplot)
library(corrr)
library(fastDummies)
library(factoextra)
library(pheatmap)
library(FactoMineR)
library(devtools)
#График ИМТ-траты в plotly (№2)
plot_ly(
data = insurance,
x = ~ bmi,
y = ~ charges,
color = ~ smoker)
plot <- insurance %>%
ggplot(aes(x = bmi, y = charges, color = smoker)) +
geom_point(size = 1.5) +
theme_light()
ggplotly(plot)
insurance_for_cor <- insurance %>%
select(is.integer | is.numeric)
insurance_cor <- cor(insurance_for_cor)
insurance_cor
## age children bmi charges
## age 1.0000000 0.04246900 0.1092719 0.29900819
## children 0.0424690 1.00000000 0.0127589 0.06799823
## bmi 0.1092719 0.01275890 1.0000000 0.19834097
## charges 0.2990082 0.06799823 0.1983410 1.00000000
corrplot(insurance_cor, method = "color", order = "alphabet", type = "upper")
corrplot.mixed(insurance_cor, lower = "color", upper = "pie", order = "AOE")
insurance_cor %>%
rplot()
# Отберем все номинативные переменные и превратим их в бинарные и удалим оригинальные колонки из датафрейма
new_insurance <- dummy_cols(insurance, select_columns = c("sex", "smoker", "region"), remove_selected_columns = TRUE)
# Стандартизуем значения переменных
insurance_scaled <- scale(new_insurance)
head(insurance_scaled)
## age bmi children charges sex_female sex_male
## [1,] -1.4382265 -0.4531506 -0.90827406 0.2984722 1.0101410 -1.0101410
## [2,] -1.5094011 0.5094306 -0.07873775 -0.9533327 -0.9892209 0.9892209
## [3,] -0.7976553 0.3831636 1.58033487 -0.7284023 -0.9892209 0.9892209
## [4,] -0.4417824 -1.3050431 -0.90827406 0.7195739 -0.9892209 0.9892209
## [5,] -0.5129570 -0.2924471 -0.90827406 -0.7765118 -0.9892209 0.9892209
## [6,] -0.5841316 -0.8073542 -0.90827406 -0.7856145 1.0101410 -1.0101410
## smoker_no smoker_yes region_northeast region_northwest region_southeast
## [1,] -1.9698501 1.9698501 -0.5650556 -0.5662062 -0.6110952
## [2,] 0.5072734 -0.5072734 -0.5650556 -0.5662062 1.6351833
## [3,] 0.5072734 -0.5072734 -0.5650556 -0.5662062 1.6351833
## [4,] 0.5072734 -0.5072734 -0.5650556 1.7648211 -0.6110952
## [5,] 0.5072734 -0.5072734 -0.5650556 1.7648211 -0.6110952
## [6,] 0.5072734 -0.5072734 -0.5650556 -0.5662062 1.6351833
## region_southwest
## [1,] 1.7648211
## [2,] -0.5662062
## [3,] -0.5662062
## [4,] -0.5662062
## [5,] -0.5662062
## [6,] -0.5662062
# И найдем дистанции
insurance_dist <- dist(insurance_scaled, method = "euclidean")
as.matrix(insurance_dist)[1:6, 1:6]
## 1 2 3 4 5 6
## 1 0.000000 5.825239 6.253322 5.747217 5.759522 4.978144
## 2 5.825239 0.000000 1.823634 4.289327 3.582563 3.361726
## 3 6.253322 1.823634 0.000000 4.663256 4.148789 3.956548
## 4 5.747217 4.289327 4.663256 0.000000 1.807952 4.583438
## 5 5.759522 3.582563 4.148789 1.807952 0.000000 4.329507
## 6 4.978144 3.361726 3.956548 4.583438 4.329507 0.000000
# Высчитываем дендрограмму
insurance_hc <- hclust(d = insurance_dist,
method = "ward.D2")
# И визуализируем
fviz_dend(insurance_hc, cex = 0.1)
fviz_dend(insurance_hc, k = 3,
cex = 0.1,
k_colors = c("#FF99FF", "#33FF99", "#00CCFF"),
color_labels_by_k = TRUE,
type = "circular"
)
clusters <- cutree(insurance_hc, k = 3)
fviz_cluster(list(data = insurance_scaled, cluster = clusters),
palette = c("#FF99FF", "#33FF99", "#00CCFF"),
ellipse.type = "convex",
repel = TRUE,
show.clust.cent = FALSE,
ggtheme = theme_light())
pheatmap(insurance_scaled)
insurance_pca <- prcomp(new_insurance, scale = T)
summary(insurance_pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.6737 1.4023 1.2417 1.1513 1.1498 1.07055 0.98583
## Proportion of Variance 0.2334 0.1639 0.1285 0.1105 0.1102 0.09551 0.08099
## Cumulative Proportion 0.2334 0.3973 0.5258 0.6363 0.7465 0.84196 0.92295
## PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.87032 0.40877 2.012e-15 1.043e-15 7.413e-16
## Proportion of Variance 0.06312 0.01392 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion 0.98608 1.00000 1.000e+00 1.000e+00 1.000e+00
Из полученных результатов видно, что 50% вариации данных объяснятся первыми тремя переменными. Посмотрим на это на графике:
fviz_eig(insurance_pca,
addlabels = T,
ylim = c(0, 30))
Посмотрим распределение переменных относительно pc1 и pc2:
fviz_pca_var(insurance_pca, col.var = "contrib")
Из графика видно, что переменные sex_male и sex_female направлены строго противоположно, что логично, так как это несовместные события. Помимо этого, образовалась группа из переменных charges и smoker_yes. Переменные регионов проживания образовали тесную группу, сосредоточенную у центра окружности. Иными словами, регион проживания вносит незначительный вклад в анализируемые главные компоненты.
Чтобы немного почистить график, оставим только топ 5 самых важных переменных:
fviz_pca_var(insurance_pca,
select.var = list(contrib = 5),
col.var = "contrib")
Чтобы понять, что стоит за pc1, pc2 и pc3, построим графики:
fviz_contrib(insurance_pca, choice = "var", axes = 1, top = 12)
fviz_contrib(insurance_pca, choice = "var", axes = 2, top = 12)
fviz_contrib(insurance_pca, choice = "var", axes = 3, top = 12)
Таким образом, переменная smoker хорошо кластериует наблюдения в оригинальном датасете. Аналогично можно судить о переменной sex.
И построим график ggbiplot, чтобы на нем одновременно изобразить переменные и наблюдения:
ggbiplot(insurance_pca,
scale = 0, alpha = 0.1) +
theme_light()
Мне захотелось визуализировать кластеризацию по полу и статусу курильщика
ggbiplot(insurance_pca,
scale=0,
groups = as.factor(insurance$sex),
ellipse = T,
alpha = 0.2) +
theme_light()
ggbiplot(insurance_pca,
scale=0,
groups = as.factor(insurance$smoker),
ellipse = T,
alpha = 0.2) +
theme_light()
# Создадим сначала переменную age_group
new_insurance_age <- new_insurance %>%
mutate(age_group = case_when(
age < 26 ~ "18-25",
age >= 26 & age < 41 ~ "26-40",
age >= 41 & age < 56 ~ "41-55",
age >= 56 ~ "56-65"
))
ggbiplot(insurance_pca,
scale=0,
groups = as.factor(new_insurance_age$age_group),
ellipse = T,
alpha = 0.2) +
theme_light()
# Создадим переменную any_child, которая будет рассказывать о том, есть ли у человека хотя бы один ребенок
new_insurance_children <- new_insurance %>%
mutate(any_child = as.character(ifelse(children >= 1, "Yes", "No")))
ggbiplot(insurance_pca,
scale=0,
groups = as.factor(new_insurance_children$any_child),
ellipse = T,
alpha = 0.2) +
theme_light()
# Создадим переменную, которая разобьет субъектов на группы в соответствии с тем, является ли их масса тела недостаточной, избыточной или нормальной
new_insurance_bmi <- new_insurance %>%
mutate(bmi_group = case_when(
bmi < 18.5 ~ "underweight",
bmi >= 18.5 & bmi < 25 ~ "normal",
bmi >= 25 ~ "overweight"
))
ggbiplot(insurance_pca,
scale=0,
groups = as.factor(new_insurance_bmi$bmi_group),
ellipse = T,
alpha = 0.2) +
theme_light()